Load data

vcs <- rio::import("data_complete_all_anonym_allZ.xlsx")
var_label(vcs$dominance) <- "Dominance"
var_label(vcs$neuro) <- "Neuroticism"
var_label(vcs$agree) <- "Agreeableness"
var_label(vcs$extra) <- "Extraversion"
var_label(vcs$openn) <- "Openness"
var_label(vcs$consc) <- "Conscientiousness"
var_label(vcs$soi_full) <- "Unrestricted sociosexuality"
var_label(vcs$f0) <- "Voice pitch"
var_label(vcs$pf) <- "Formants"
vcs$sex_c <- vcs$sex
vcs$sex <- factor(if_else(vcs$sex == 1, "male", "female"))
contrasts(vcs$sex) <- contr.helmert(2)
var_label(vcs$sex) <- "Sex"
set.seed(1)
var_label(vcs$age) <- "Age"
vcs <- vcs %>% 
  mutate(
    age = if_else(dataset == 9, 20, age)/10,
    age_se = if_else(dataset == 9, 3, 0.5)/10
  )
warmup <- 2000
iter <- warmup + 2000
chains <- 4
control <- list(adapt_delta = 0.99)
priors <- c(
  prior(normal(0, 3), class = b)
)

variable_labels <- c("Intercept"= "Intercept",  "sex1" = "Sex [male]", "bsp_meageage_se" = "Age±SE", "f0" = "Voice pitch (f0)", "pf" = "Formants (f1-f4)", "age" = "Age")
effect_labels <- c("b_Intercept"= "Intercept",  "b_sex1" = "Sex [male]", "bsp_meageage_se" = "Age±SE", "b_f0" = "Voice pitch (f0)", "b_pf" = "Formants (f1-f4)", "b_age" = "Age")

Missingness patterns

codebook::md_pattern(vcs %>% select(dominance, sex, age, f0, pf, dataset))
## # A tibble: 4 x 7
##   description                       age    pf    f0 dominance var_miss n_miss
##   <chr>                           <dbl> <dbl> <dbl>     <dbl>    <dbl>  <dbl>
## 1 Missing values per variable         4     4     7      1256     1271   1271
## 2 Missing values in 1 variables       1     1     1         0        1   1245
## 3 Missing values in 0 variables       1     1     1         1        0    985
## 4 3 other, less frequent patterns     2     2     1         0        7     11
vcs <- vcs %>% filter(!is.na(f0), !is.na(pf), !is.na(age))

Hypothesis 1:

Participants with lower voice pitch will have a more dominant personality.

Interpretation: The data are are consistent with substantial linear negative relationships between f0 and dominance. The 89% HDI for f0 falls entirely outside the ROPE, but the HDI for pf falls almost entirely within it. This evidence is consistent with a non-negligible association, where people with deeper voices are more dominant (after adjusting for age and gender). Further research is needed to verify if the association with pf is truly negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).

## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing missing values (geom_point).

ggplot(vcs, aes(pf, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing missing values (geom_point).

ggplot(vcs, aes(age*10, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing missing values (geom_point).

Decisions upon visual diagnosis: Visual diagnostic show no clear sign of nonlinearities. We will fit linear effects for f0, pf, and age. It is not necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is not included in these analyses. Visually, there could be an interaction between sex and f0.

Models

h1_simple <- brm(dominance ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/dominance/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/dominance/h1_simple.rds'
## Automatically saving the model object in 'models/dominance/h1_simple.rds'
## Automatically saving the model object in 'models/dominance/h1_simple.rds'
h1_after_diagnosis <- brm(dominance ~ f0*sex + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, cores = chains, chains = chains,
          prior = priors, 
          control = control, save_mevars = TRUE,
          file = "models/dominance/h1_nonlinear") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/dominance/h1_nonlinear.rds'
## Automatically saving the model object in 'models/dominance/h1_nonlinear.rds'
## Automatically saving the model object in 'models/dominance/h1_nonlinear.rds'
# conditional_smooths(h1_after_diagnosis)

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_simple           0.0       0.0   
## h1_after_diagnosis -1.0       0.2
h1_sex_mod_dataset_varying <- brm(dominance ~ f0 + sex + age + 
                                    (1 + f0 || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/dominance/h1_sex_mod_dataset_varying") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/dominance/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/dominance/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/dominance/h1_sex_mod_dataset_varying.rds'
h1_after_diagnosis_melsm <- brm(
  bf(dominance ~ f0 + sex + age + (1 | dataset),
     sigma ~ f0 + sex + age + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/dominance/h1_after_diagnosis_melsm") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/dominance/h1_after_diagnosis_melsm.rds'
## Automatically saving the model object in 'models/dominance/h1_after_diagnosis_melsm.rds'
## Automatically saving the model object in 'models/dominance/h1_after_diagnosis_melsm.rds'
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm)
compare_models
##                            elpd_diff se_diff
## h1_after_diagnosis_melsm    0.0       0.0   
## h1_simple                  -0.2       4.0   
## h1_after_diagnosis         -1.2       4.0   
## h1_sex_mod_dataset_varying -2.9       4.4
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  Dominance
Predictors Estimates CI (95%)
Intercept -0.08 -0.80;0.60
Voice pitch (f0) -0.27 -0.46;-0.08
Formants (f1-f4) -0.02 -0.20;0.16
Sex [male] -0.31 -0.52;-0.08
Age±SE 0.00 -0.01;0.02
Random Effects
σ2 0.05
Ï„00 0.96
ICC 0.05
N dataset 4
Observations 985
Marginal R2 / Conditional R2 0.069 / 0.069 
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
## Possible multicollinearity between b_sex1 and b_f0 (r = 0.76). This might lead to inappropriate results. See 'Details' in '?equivalence_test'.
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Undecided 32% [-0.60; 0.41]
Voice pitch (f0) Rejected 0% [-0.42;-0.11]
Formants (f1-f4) Undecided 81% [-0.16; 0.13]
Sex [male] Rejected 0% [-0.50;-0.14]
Age±SE Accepted 100% [-0.01; 0.01] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Picking joint bandwidth of 0.00905
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details

summary(h1_simple)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 5 divergent transitions after warmup. Increasing adapt_delta
## above 0.99 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.52      0.43     0.16     1.66 1.00     1854     2414
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      -0.08      0.37    -0.80     0.60 1.00     1470     1979
## f0             -0.27      0.10    -0.46    -0.08 1.00     4886     5623
## pf             -0.02      0.09    -0.20     0.16 1.00     7512     6543
## sex1           -0.31      0.11    -0.52    -0.08 1.00     4367     5610
## meageage_se     0.00      0.01    -0.01     0.02 1.00     7099     5994
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.93     1.01 1.00     8321     5117
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta
## above 0.99 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 * sex + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.53      0.43     0.16     1.63 1.00     2619     3493
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      -0.05      0.37    -0.79     0.67 1.00     3111     3051
## f0             -0.26      0.10    -0.45    -0.07 1.00     6625     5759
## sex1           -0.30      0.12    -0.52    -0.07 1.00     5716     5773
## pf             -0.02      0.09    -0.19     0.16 1.00     9331     6207
## f0:sex1         0.02      0.10    -0.17     0.21 1.00     8589     6232
## meageage_se     0.00      0.01    -0.01     0.02 1.00     9747     6170
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.93     1.02 1.00    10519     5550
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_after_diagnosis_melsm)

# conditional_smooths(h1_after_diagnosis_melsm)
summary(h1_sex_mod_dataset_varying)
## Warning: There were 4 divergent transitions after warmup. Increasing adapt_delta
## above 0.99 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 + pf + f0 * sex + me(age, age_se) + (1 + f0 * sex || dataset) 
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.50      0.47     0.03     1.78 1.00     3660     3397
## sd(f0)            0.26      0.35     0.01     1.20 1.00     3195     3060
## sd(sex1)          0.34      0.43     0.01     1.51 1.00     2747     3104
## sd(f0:sex1)       0.37      0.39     0.02     1.41 1.00     2894     3588
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      -0.09      0.39    -0.84     0.68 1.00     5857     4642
## f0             -0.25      0.26    -0.71     0.20 1.00     3824     2480
## pf             -0.03      0.09    -0.21     0.16 1.00    11668     8835
## sex1           -0.29      0.30    -0.95     0.29 1.00     3959     2492
## f0:sex1        -0.00      0.28    -0.56     0.54 1.00     4020     3426
## meageage_se     0.00      0.01    -0.01     0.02 1.00    11857     9601
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.93     1.01 1.00    13009     8654
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis_melsm)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta
## above 0.99 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: dominance ~ f0 * sex + pf + sex + me(age, age_se) + (1 | dataset) 
##          sigma ~ f0 * sex + pf + sex + me(age, age_se) + (1 | dataset)
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.52      0.39     0.17     1.58 1.00     3584     4787
## sd(sigma_Intercept)     0.15      0.17     0.01     0.62 1.00     2706     3542
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -0.07      0.35    -0.79     0.60 1.00     3595     3307
## sigma_Intercept       0.04      0.19    -0.34     0.38 1.00     4339     5652
## f0                   -0.29      0.10    -0.48    -0.09 1.00     7558     8924
## sex1                 -0.30      0.12    -0.54    -0.07 1.00     6195     7846
## pf                    0.00      0.09    -0.16     0.18 1.00    11479     9318
## f0:sex1               0.04      0.10    -0.15     0.23 1.00    12251     9131
## sigma_f0              0.10      0.07    -0.04     0.25 1.00     7277     8700
## sigma_sex1            0.13      0.09    -0.04     0.30 1.00     6632     8402
## sigma_pf              0.03      0.07    -0.11     0.16 1.00    11282     9463
## sigma_f0:sex1        -0.01      0.07    -0.15     0.13 1.00    10718     8125
## meageage_se           0.00      0.01    -0.01     0.02 1.00    11342     9225
## sigma_meageage_se    -0.00      0.01    -0.01     0.01 1.00     4908     7372
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_simple 0.0699913 0.0147882 0.0427074 0.1001835 0.0529997
h1_after_diagnosis 0.0705955 0.0149998 0.0434176 0.1014367 0.0511312
h1_after_diagnosis_melsm 0.0723500 0.0145552 0.0455458 0.1020524 0.0509098
h1_sex_mod_dataset_varying 0.0754052 0.0149069 0.0474739 0.1061981 0.0476857

Hypothesis 2:

Participants with lower voice pitch will score higher on agreeableness.

Interpretation: The credible intervals for f0 and pf include zero, but the 89% HDI does not fall entirely within the ROPE. This evidence is consistent with a negligible association between f0/pf and agreeableness. Further research is needed to verify if the associations are truly negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour =  "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 796 rows containing non-finite values (stat_smooth).

## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).

ggplot(vcs, aes(pf, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).

ggplot(vcs, aes(age*10, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).

Decisions upon visual diagnosis: Visual diagnostic shows signs of nonlinearity for f0 and pf, but no clear interactions by sex. We will fit splines for f0, and pf. It is necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is included in these analyses.

Models

h1_simple <- brm(agree ~ f0 + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/agree/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/agree/h1_simple.rds'
## Automatically saving the model object in 'models/agree/h1_simple.rds'
## Automatically saving the model object in 'models/agree/h1_simple.rds'
h1_after_diagnosis <- brm(agree ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, cores = chains, chains = chains,
          prior = priors, 
          control = control, save_mevars = TRUE,
          file = "models/agree/h1_nonlinear") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/agree/h1_nonlinear.rds'
## Automatically saving the model object in 'models/agree/h1_nonlinear.rds'
## Automatically saving the model object in 'models/agree/h1_nonlinear.rds'
conditional_smooths(h1_after_diagnosis)

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_after_diagnosis  0.0       0.0   
## h1_simple          -3.8       3.1
h1_sex_mod_dataset_varying <- brm(agree ~ f0 + pf + f0 * sex + me(age, age_se) + 
                                    (1 + f0 + pf + sex || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/agree/h1_sex_mod_dataset_varying") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/agree/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/agree/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/agree/h1_sex_mod_dataset_varying.rds'
h1_after_diagnosis_melsm <- brm(
  bf(agree ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset),
     sigma ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/agree/h1_after_diagnosis_melsm") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/agree/h1_after_diagnosis_melsm.rds'
## Automatically saving the model object in 'models/agree/h1_after_diagnosis_melsm.rds'
## Automatically saving the model object in 'models/agree/h1_after_diagnosis_melsm.rds'
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm)
compare_models
##                            elpd_diff se_diff
## h1_after_diagnosis_melsm     0.0       0.0  
## h1_after_diagnosis          -7.0       5.1  
## h1_sex_mod_dataset_varying  -7.5       6.2  
## h1_simple                  -10.9       5.9
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  Agreeableness
Predictors Estimates CI (95%)
Intercept 0.09 -0.37;0.59
Voice pitch (f0) 0.01 -0.13;0.16
Formants (f1-f4) 0.09 -0.06;0.24
Sex [male] 0.02 -0.15;0.20
Age±SE -0.01 -0.02;0.00
Random Effects
σ2 0.11
Ï„00 0.90
ICC 0.11
N dataset 7
Observations 1434
Marginal R2 / Conditional R2 0.117 / 0.117 
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
## Possible multicollinearity between b_sex1 and b_f0 (r = 0.73). This might lead to inappropriate results. See 'Details' in '?equivalence_test'.
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Undecided 35% [-0.28; 0.47]
Voice pitch (f0) Undecided 91% [-0.11; 0.13]
Formants (f1-f4) Undecided 58% [-0.04; 0.20]
Sex [male] Undecided 80% [-0.13; 0.16]
Age±SE Accepted 100% [-0.01; 0.00] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Picking joint bandwidth of 0.00726
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details

summary(h1_simple)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ f0 + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.48      0.21     0.25     0.99 1.00     1992     3446
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.09      0.25    -0.37     0.59 1.00     2017     3031
## f0              0.01      0.07    -0.13     0.16 1.00     6296     5690
## pf              0.09      0.08    -0.06     0.24 1.00     7233     5707
## sex1            0.02      0.09    -0.15     0.20 1.00     5467     5890
## meageage_se    -0.01      0.01    -0.02     0.00 1.00     8114     6260
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.94      0.02     0.91     0.98 1.00    11815     6095
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Smooth Terms: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sf0_1)     1.05      0.65     0.23     2.73 1.00     3297     3480
## sds(spf_1)     0.65      0.50     0.03     1.91 1.00     3618     3082
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.48      0.20     0.24     1.01 1.00     2375     2765
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.11      0.24    -0.35     0.60 1.00     2712     3234
## sex1            0.13      0.13    -0.11     0.39 1.00     7526     6076
## sf0_1          -0.24      1.58    -3.54     2.84 1.00     6436     5755
## spf_1           0.60      1.26    -2.04     3.33 1.00     5442     4987
## meageage_se    -0.01      0.01    -0.02     0.00 1.00    10169     6084
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.94      0.02     0.90     0.97 1.00    14100     6179
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_after_diagnosis_melsm)

# conditional_smooths(h1_after_diagnosis_melsm)
summary(h1_sex_mod_dataset_varying)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ f0 + pf + f0 * sex + me(age, age_se) + (1 + f0 + pf + sex || dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.53      0.23     0.26     1.14 1.00     3925     5371
## sd(f0)            0.15      0.12     0.01     0.44 1.00     3694     4711
## sd(pf)            0.14      0.12     0.01     0.44 1.00     4726     5476
## sd(sex1)          0.15      0.16     0.01     0.57 1.00     3137     4980
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.22      0.28    -0.34     0.77 1.00     3801     5048
## f0              0.04      0.11    -0.17     0.25 1.00     6903     5540
## pf              0.04      0.11    -0.18     0.24 1.00     7535     6008
## sex1            0.07      0.15    -0.23     0.35 1.00     5261     4665
## f0:sex1         0.13      0.08    -0.04     0.30 1.00    10675     8696
## meageage_se    -0.01      0.01    -0.02     0.01 1.00    11214     8721
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.94      0.02     0.90     0.97 1.00    16491     8626
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis_melsm)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: agree ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset) 
##          sigma ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset)
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Smooth Terms: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sf0_1)           1.24      0.77     0.27     3.20 1.00     3949     5731
## sds(spf_1)           0.60      0.49     0.03     1.88 1.00     5289     6287
## sds(sigma_sf0_1)     0.49      0.49     0.01     1.79 1.00     3876     6866
## sds(sigma_spf_1)     0.77      0.61     0.03     2.29 1.00     3279     5300
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.47      0.19     0.24     0.95 1.00     4396     6685
## sd(sigma_Intercept)     0.12      0.07     0.02     0.28 1.00     3323     3103
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.13      0.23    -0.32     0.58 1.00     4551     5867
## sigma_Intercept       0.18      0.12    -0.05     0.41 1.00     7990     7681
## sex1                  0.16      0.13    -0.08     0.43 1.00    10684     9656
## sigma_sex1           -0.08      0.09    -0.29     0.09 1.00     9147     7879
## sf0_1                -0.28      1.70    -3.76     3.02 1.00     9243     8390
## spf_1                 0.71      1.21    -1.85     3.15 1.00     8565     6811
## sigma_sf0_1          -0.28      1.31    -2.61     3.01 1.00     6575     4527
## sigma_spf_1          -1.32      1.73    -5.48     1.38 1.00     5933     6055
## meageage_se          -0.01      0.01    -0.02     0.00 1.00    15953     8469
## sigma_meageage_se    -0.01      0.00    -0.02    -0.00 1.00     9315     8745
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_after_diagnosis_melsm 0.1264509 0.0160206 0.0953042 0.1580384 0.1098105
h1_after_diagnosis 0.1266616 0.0157463 0.0965032 0.1587202 0.1094989
h1_sex_mod_dataset_varying 0.1301130 0.0154954 0.1003210 0.1609634 0.1085954
h1_simple 0.1176955 0.0153079 0.0881099 0.1484517 0.1045678

Hypothesis 3:

Participants with lower voice pitch will score lower on neuroticism.

Interpretation: The data are are consistent with linear positive relationship between f0 and neuroticism. The 89% HDI for f0 excludes zero but partly inside the ROPE. 69% of the HDI for pf falls within the ROPE. This evidence is consistent with an association, where people with deeper voices have lower neuroticism (after adjusting for age and gender), but further research is needed to test whether the associations with f0 and pf are negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, neuro)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour =  "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 796 rows containing non-finite values (stat_smooth).

## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).

ggplot(vcs, aes(pf, neuro)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).

ggplot(vcs, aes(age*10, neuro)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).

Decisions upon visual diagnosis: Visual diagnostic shows signs of a potential interaction by sex f0 and a potential interaction by sex coupled with nonlinearity for pf. We will fit an interaction for f0, and separate splines by sex for pf. It is necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is included in these analyses.

Models

h1_simple <- brm(neuro ~ f0 + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/neuro/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/neuro/h1_simple.rds'
## Automatically saving the model object in 'models/neuro/h1_simple.rds'
## Automatically saving the model object in 'models/neuro/h1_simple.rds'
h1_after_diagnosis <- brm(neuro ~ f0 * sex + s(pf, by = sex) + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, cores = chains, chains = chains,
          prior = priors, 
          control = control, save_mevars = TRUE,
          file = "models/neuro/h1_nonlinear") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/neuro/h1_nonlinear.rds'
## Automatically saving the model object in 'models/neuro/h1_nonlinear.rds'
## Automatically saving the model object in 'models/neuro/h1_nonlinear.rds'
conditional_smooths(h1_after_diagnosis)

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_after_diagnosis  0.0       0.0   
## h1_simple          -0.6       2.4
h1_sex_mod_dataset_varying <- brm(neuro ~ f0 + sex + me(age, age_se) + 
                                    (1 + f0 + sex + age || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/neuro/h1_sex_mod_dataset_varying") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/neuro/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/neuro/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/neuro/h1_sex_mod_dataset_varying.rds'
# agree <- brm(
#   bf(neuro ~ f0 + sex + me(age, age_se) + (1 | dataset),
#      sigma ~ f0 + sex + me(age, age_se) + (1 | dataset)), data = vcs,
#           prior = priors,
#           iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
#           control = control, save_mevars = TRUE,
#           file = "models/neuro/agree") %>% 
#   add_criterion("loo") %>% 
#   add_criterion("bayes_R2") %>% 
#   add_criterion("loo_R2") 
#   

compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying
                      # ,h1_after_diagnosis_melsm
                      )
compare_models
##                            elpd_diff se_diff
## h1_sex_mod_dataset_varying  0.0       0.0   
## h1_after_diagnosis         -2.1       3.3   
## h1_simple                  -2.7       4.0
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying
                      # ,h1_after_diagnosis_melsm=h1_after_diagnosis_melsm
            ) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  Neuroticism
Predictors Estimates CI (95%)
Intercept 0.04 -0.28;0.35
Voice pitch (f0) 0.16 0.01;0.30
Formants (f1-f4) -0.07 -0.21;0.06
Sex [male] -0.15 -0.33;0.02
Age±SE -0.02 -0.14;0.09
Random Effects
σ2 0.01
Ï„00 0.99
ICC 0.01
N dataset 7
Observations 1434
Marginal R2 / Conditional R2 0.088 / 0.088 
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
## Possible multicollinearity between b_sex1 and b_f0 (r = 0.77). This might lead to inappropriate results. See 'Details' in '?equivalence_test'.
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Undecided 51% [-0.22; 0.30]
Voice pitch (f0) Undecided 19% [ 0.03; 0.27]
Formants (f1-f4) Undecided 69% [-0.18; 0.04]
Sex [male] Undecided 25% [-0.29;-0.01]
Age±SE Undecided 97% [-0.11; 0.07] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Picking joint bandwidth of 0.00851
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details

summary(h1_simple)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: neuro ~ f0 + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.14      0.08     0.03     0.34 1.00     1546     1455
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.04      0.16    -0.28     0.35 1.00     2250     3563
## f0              0.16      0.07     0.01     0.30 1.00     4243     5581
## pf             -0.07      0.07    -0.21     0.06 1.00     4646     5858
## sex1           -0.15      0.09    -0.33     0.02 1.00     3609     4796
## meageage_se    -0.02      0.06    -0.14     0.09 1.00     2512     3192
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.96      0.02     0.92     0.99 1.00     9149     5924
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: neuro ~ f0 * sex + s(pf, by = sex) + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Smooth Terms: 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(spfsexfemale_1)     0.63      0.59     0.02     2.16 1.00     4619     3779
## sds(spfsexmale_1)       1.22      0.94     0.07     3.67 1.00     3335     3158
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.15      0.08     0.04     0.35 1.00     3066     3700
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          -0.22      0.21    -0.64     0.19 1.00     4524     5308
## f0                  0.12      0.08    -0.04     0.27 1.00     7955     6568
## sex1               -0.26      0.13    -0.51    -0.02 1.00     5197     5527
## f0:sex1            -0.16      0.08    -0.31    -0.01 1.00    10869     6230
## spf:sexfemale_1    -0.24      1.50    -3.27     3.03 1.00     4854     4458
## spf:sexmale_1      -0.13      2.07    -4.28     4.26 1.00     6130     5760
## meageage_se        -0.01      0.06    -0.12     0.10 1.00     5308     5097
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.95      0.02     0.92     0.99 1.00    10209     5305
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_after_diagnosis_melsm)

# conditional_smooths(h1_after_diagnosis_melsm)
summary(h1_sex_mod_dataset_varying)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: neuro ~ f0 + pf * sex + f0 * sex + me(age, age_se) + (1 + f0 * sex + pf * sex + sex || dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.24      0.17     0.02     0.67 1.00     3540     4464
## sd(f0)            0.13      0.12     0.00     0.43 1.00     4402     6207
## sd(sex1)          0.16      0.14     0.01     0.52 1.00     4228     5443
## sd(pf)            0.19      0.17     0.01     0.64 1.00     4928     6626
## sd(f0:sex1)       0.14      0.12     0.01     0.45 1.00     5134     6440
## sd(sex1:pf)       0.31      0.24     0.02     0.89 1.00     3015     4627
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      -0.21      0.22    -0.65     0.21 1.00     9960     8724
## f0              0.13      0.11    -0.08     0.35 1.00     8530     6968
## pf             -0.12      0.15    -0.42     0.18 1.00     7696     8006
## sex1           -0.20      0.14    -0.47     0.08 1.00     8816     7551
## pf:sex1        -0.10      0.19    -0.47     0.26 1.00     5858     4908
## f0:sex1        -0.14      0.11    -0.36     0.07 1.00     9590     7033
## meageage_se     0.01      0.06    -0.11     0.12 1.00    11138     8693
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.95      0.02     0.92     0.99 1.00    19230     8270
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# summary(h1_after_diagnosis_melsm)
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_sex_mod_dataset_varying 0.1065406 0.0144428 0.0791614 0.1352072 0.0785213
h1_after_diagnosis 0.0943398 0.0139157 0.0682566 0.1215960 0.0762639
h1_simple 0.0880062 0.0134624 0.0622406 0.1155778 0.0753931

Hypothesis 4:

Participants with lower voice pitch will report having a more unrestricted sociosexuality.

Interpretation: The data are are consistent with substantial linear negative relationship between f0 and unrestricted sociosexuality. The 89% HDI for f0 falls entirely outside the ROPE, but the HDI for pf falls almost entirely within it. This evidence is consistent with a non-negligible association, where people with deeper voices have a less restricted sociosexuality (after adjusting for age and gender). Further research is needed to verify if the association with pf is truly negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, soi_full)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 235 rows containing non-finite values (stat_smooth).

## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing missing values (geom_point).

ggplot(vcs, aes(pf, soi_full)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing missing values (geom_point).

ggplot(vcs, aes(age*10, soi_full)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing missing values (geom_point).

Decisions upon visual diagnosis: We will fit splines for f0, pf, and age. It seems likely that the nonlinearities in f0 and pf will not both be apparent after adjusting for one another. There is no evidence in the plots that there is an interaction with sex, except for age.

Models

h1_simple <- brm(soi_full ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/soi_full/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/soi_full/h1_simple.rds'
## Automatically saving the model object in 'models/soi_full/h1_simple.rds'
## Automatically saving the model object in 'models/soi_full/h1_simple.rds'
h1_after_diagnosis <- brm(soi_full ~  s(f0) + s(pf) + sex + s(age, by = sex) + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, cores = chains, chains = chains,
          prior = priors, 
          control = control, save_mevars = TRUE,
          file = "models/soi_full/h1_nonlinear") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/soi_full/h1_nonlinear.rds'
## Automatically saving the model object in 'models/soi_full/h1_nonlinear.rds'
## Automatically saving the model object in 'models/soi_full/h1_nonlinear.rds'
conditional_smooths(h1_after_diagnosis)

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_after_diagnosis  0.0       0.0   
## h1_simple          -3.5       3.6
h1_sex_mod_dataset_varying <- brm(soi_full ~ f0 + sex + s(age, by = sex) + 
                                    (1 + f0 + sex + age || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/soi_full/h1_sex_mod_dataset_varying") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Automatically saving the model object in 'models/soi_full/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/soi_full/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/soi_full/h1_sex_mod_dataset_varying.rds'
h1_after_diagnosis_melsm <- brm(
  bf(soi_full ~  f0 + sex + s(age, by = sex) + (1 | dataset),
     sigma ~ f0 + sex + s(age, by = sex) + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/soi_full/h1_after_diagnosis_melsm") %>%
  add_criterion("loo") %>%
  add_criterion("bayes_R2") %>%
  add_criterion("loo_R2")
## Automatically saving the model object in 'models/soi_full/h1_after_diagnosis_melsm.rds'
## Automatically saving the model object in 'models/soi_full/h1_after_diagnosis_melsm.rds'
## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
## Automatically saving the model object in 'models/soi_full/h1_after_diagnosis_melsm.rds'
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying, h1_after_diagnosis_melsm)
compare_models
##                            elpd_diff se_diff
## h1_after_diagnosis_melsm     0.0       0.0  
## h1_after_diagnosis          -8.2       5.5  
## h1_sex_mod_dataset_varying -11.4       5.4  
## h1_simple                  -11.7       6.0
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying, h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  Unrestricted
sociosexuality
Predictors Estimates CI (95%)
Intercept -0.20 -0.60;0.19
Voice pitch (f0) -0.29 -0.41;-0.17
Formants (f1-f4) -0.01 -0.14;0.12
Sex [male] -0.03 -0.17;0.12
Age 0.09 -0.01;0.18
Random Effects
σ2 0.08
Ï„00 0.93
ICC 0.07
N dataset 9
Observations 1995
Marginal R2 / Conditional R2 0.160 / 0.160 
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Undecided 26% [-0.51; 0.13]
Voice pitch (f0) Rejected 0% [-0.39;-0.19]
Formants (f1-f4) Undecided 96% [-0.12; 0.09]
Sex [male] Undecided 89% [-0.14; 0.09]
Age Undecided 60% [ 0.02; 0.16] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Picking joint bandwidth of 0.0073
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details

summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: soi_full ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1995) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.46      0.15     0.27     0.85 1.00     1910     3239
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.20      0.20    -0.60     0.19 1.00     2156     3365
## f0           -0.29      0.06    -0.41    -0.17 1.00     4808     4766
## pf           -0.01      0.06    -0.14     0.12 1.00     6259     5195
## sex1         -0.03      0.07    -0.17     0.12 1.00     4258     4948
## age           0.09      0.05    -0.01     0.18 1.00     7240     5169
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.92      0.01     0.89     0.95 1.00     7124     4467
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: soi_full ~ s(f0) + s(pf) + sex + s(age, by = sex) + (1 | dataset) 
##    Data: vcs (Number of observations: 1995) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Smooth Terms: 
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sds(sf0_1)               1.62      1.08     0.08     4.11 1.00     1484
## sds(spf_1)               0.67      0.58     0.02     2.15 1.00     3178
## sds(sagesexfemale_1)     1.06      0.68     0.19     2.80 1.00     2965
## sds(sagesexmale_1)       0.49      0.46     0.02     1.72 1.00     3657
##                      Tail_ESS
## sds(sf0_1)               2166
## sds(spf_1)               4334
## sds(sagesexfemale_1)     2570
## sds(sagesexmale_1)       4202
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.44      0.14     0.25     0.80 1.00     2878     4120
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            0.02      0.16    -0.31     0.33 1.00     1905     3744
## sex1                 0.02      0.12    -0.20     0.26 1.00     6673     6118
## sf0_1               -1.68      2.06    -5.25     3.00 1.00     4230     5063
## spf_1                0.02      1.26    -2.52     2.75 1.00     6463     5365
## sage:sexfemale_1     0.88      1.70    -2.67     4.38 1.00     5886     5216
## sage:sexmale_1       0.40      1.05    -1.82     2.67 1.00     5799     4865
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.92      0.01     0.89     0.95 1.00    12164     5749
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_simple)

conditional_smooths(h1_after_diagnosis)

summary(h1_sex_mod_dataset_varying)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: soi_full ~ f0 + pf * sex + s(age, by = sex) + (1 + pf + f0 + sex || dataset) 
##    Data: vcs (Number of observations: 1995) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Smooth Terms: 
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sds(sagesexfemale_1)     1.01      0.66     0.17     2.68 1.00     5601
## sds(sagesexmale_1)       0.62      0.55     0.03     2.04 1.00     4981
##                      Tail_ESS
## sds(sagesexfemale_1)     4654
## sds(sagesexmale_1)       6139
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.45      0.16     0.24     0.84 1.00     5523     7210
## sd(pf)            0.13      0.11     0.01     0.39 1.00     4931     6390
## sd(f0)            0.10      0.08     0.00     0.30 1.00     4961     6208
## sd(sex1)          0.13      0.13     0.00     0.47 1.00     4191     5638
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -0.03      0.17    -0.39     0.31 1.00     4576     6173
## f0                  -0.29      0.08    -0.44    -0.14 1.00    10385     7627
## pf                  -0.02      0.10    -0.22     0.18 1.00     9863     8486
## sex1                -0.02      0.12    -0.26     0.21 1.00     8065     7178
## pf:sex1             -0.03      0.09    -0.21     0.14 1.00    11580     8726
## sage:sexfemale_1     0.70      1.67    -2.84     4.02 1.00     8781     8376
## sage:sexmale_1       0.62      1.22    -1.83     3.28 1.00     7634     6686
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.92      0.01     0.89     0.95 1.00    18769     8603
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis_melsm)
## Warning: There were 10 divergent transitions after warmup. Increasing
## adapt_delta above 0.99 may help. See http://mc-stan.org/misc/
## warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: soi_full ~ f0 + pf + sex + s(age, by = sex) + (1 | dataset) 
##          sigma ~ f0 + pf + sex + s(age, by = sex) + (1 | dataset)
##    Data: vcs (Number of observations: 1995) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Smooth Terms: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sds(sagesexfemale_1)           1.11      0.67     0.26     2.79 1.00     6683
## sds(sagesexmale_1)             0.44      0.42     0.01     1.58 1.00     6492
## sds(sigma_sagesexfemale_1)     0.63      0.61     0.02     2.25 1.00     3904
## sds(sigma_sagesexmale_1)       0.40      0.43     0.01     1.59 1.00     5024
##                            Tail_ESS
## sds(sagesexfemale_1)           6966
## sds(sagesexmale_1)             7184
## sds(sigma_sagesexfemale_1)     7777
## sds(sigma_sagesexmale_1)       6947
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.45      0.15     0.26     0.82 1.00     4973     7393
## sd(sigma_Intercept)     0.17      0.07     0.08     0.33 1.00     3835     6626
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                  0.02      0.16    -0.29     0.35 1.00     3550
## sigma_Intercept           -0.12      0.06    -0.25     0.01 1.00     4597
## f0                        -0.29      0.06    -0.41    -0.17 1.00    16252
## pf                         0.01      0.07    -0.12     0.14 1.00    18032
## sex1                      -0.00      0.07    -0.15     0.14 1.00    12734
## sigma_f0                  -0.00      0.05    -0.09     0.09 1.00    16967
## sigma_pf                   0.03      0.05    -0.08     0.13 1.00    17853
## sigma_sex1                 0.05      0.06    -0.06     0.16 1.00    12697
## sage:sexfemale_1           0.80      1.74    -2.91     4.29 1.00    10802
## sage:sexmale_1             0.25      0.98    -1.86     2.29 1.00     9187
## sigma_sage:sexfemale_1    -0.46      1.75    -4.77     2.73 1.00     7930
## sigma_sage:sexmale_1      -0.11      1.03    -2.06     2.32 1.00     7920
##                        Tail_ESS
## Intercept                  4945
## sigma_Intercept            5941
## f0                         8933
## pf                         9395
## sex1                       8948
## sigma_f0                   9549
## sigma_pf                   9052
## sigma_sex1                 8763
## sage:sexfemale_1           8699
## sage:sexmale_1             7383
## sigma_sage:sexfemale_1     5720
## sigma_sage:sexmale_1       5717
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_after_diagnosis 0.1688764 0.0139599 0.1416838 0.1963915 0.1537734
h1_after_diagnosis_melsm 0.1679133 0.0121578 0.1442473 0.1917065 0.1528378
h1_sex_mod_dataset_varying 0.1688764 0.0140204 0.1418163 0.1967623 0.1512160
h1_simple 0.1598969 0.0137611 0.1322719 0.1866512 0.1509295